#ifndef DYN_ODE_DOPRI853_H
#define DYN_ODE_DOPRI853_H

#include "dyn_ode.h"

/*! \brief Implementation of the Dormand-Prince
 * method of order 8(5-3).
 *
 * This method uses a eighth-order Runge-Kutta
 * order along with embedded Runge-Kutta methods
 * for error estimation. We try to minimize this
 * error with a PI stepsize control.
 */


class DOPRI853 : public ODE
{
public:
    /*! The constructor setups the integrator. */
    DOPRI853(MultiFunctor& func,
             colvec _initCond,
             double _start,
             double _end,
             double _initStepsize,
             double _minStepsize,
             double _maxStepsize,
             double _absTol,
             double _relTol,
             int _maxIterations);

    /*! @name Inherited Virtual Functions.
     * Implementation of the inherited pure
     * virtual functions.
     */
    //@{
    int integrate();
    colvec next(colvec previousStep);
    //@}

    /*! @name Accessor Functions
     * We access the old time step and
     * the PI control scheme variables.
     */
    //@{
    double getOldStep(){return oldStep;}
    double getBeta(){return beta;}
    double getAlpha(){return alpha;}

    void setOldStep(double _oldStep){oldStep=_oldStep;}
    void setBeta(double _beta){beta=_beta;}
    void setAlaph(double _alpha){alpha=_alpha;}
    //@}

protected:
    /*! Size of the step we have
     * previously taken.
     */
    double oldStep;

    /*! @name PI Stepsize Control.
     * We hold the last error we have in memory
     * to make the stepsize correction.
     */
    //@{
    double errorOld=1e-4;
    double beta=0.05;
    double alpha=0.125-0.75*beta;
    //@}

    /*! @name Butcher Tableau
     * The following constants form
     * the Butcher tableau of the method.
     */
    //@{
    /*! @name Runge-Kutta matrix.*/
    //@{
    /*! Eases up expressions. */
    const double sq6 = sqrt(6.);

    const double a21 = (12.-2.*sq6)/135.;
    const double a31 = (6.-sq6)/180.;
    const double a32 = (6.-sq6)/60.;
    const double a41 = (6.-sq6)/120.;
    const double a43 = (6.-sq6)/40.;
    const double a51 = (462.+107.*sq6)/3000.;
    const double a53 = -(402.+197.*sq6)/1000.;
    const double a54 = (168.+73.*sq6)/375.;
    const double a61 = 1./27.;
    const double a64 = (16.+sq6)/108.;
    const double a65 = (16.-sq6)/108.;
    const double a71 = 19./512.;
    const double a74 = (118.+23.*sq6)/1024.;
    const double a75 = (118.-23.*sq6)/1024.;
    const double a76 = -9./512.;
    const double a81 = 13772./371293.;
    const double a84 = (51544.+4784.*sq6)/371293.;
    const double a85 = (51544.-4784.*sq6)/371293.;
    const double a86 = -5688./371293.;
    const double a87 = 3072./371293.;
    const double a91 = 58656157643./93983540625.;
    const double a94 = -(1324889724104.+318801444819.*sq6)/626556937500.;
    const double a95 = -(1324889724104.-318801444819.*sq6)/626556937500.;
    const double a96 = 96044563816./3480871875.;
    const double a97 = 5682451879168./281950621875.;
    const double a98 = -165125654./3796875.;
    const double a101 = 8909899./18653125.;
    const double a104 = -(4521408.+1137963.*sq6)/2937500.;
    const double a105 = -(4521408.-1137963.*sq6)/2937500.;
    const double a106 = 96663078./4553125.;
    const double a107 = 2107245056./137915625.;
    const double a108 = -4913652016./147609375.;
    const double a109 = -78894270./3880452869.;
    const double a111 = -20401265806./21769653311.;
    const double a114 = (354216.+94326.*sq6)/112847.;
    const double a115 = (354216.-94326.*sq6)/112847.;
    const double a116 = -43306765128./5313852383.;
    const double a117 = -20866708358144./1126708119789.;
    const double a118 = 14886003438020./654632330667.;
    const double a119 = 35290686222309375./14152473387134411.;
    const double a1110 = -1477884375./485066827.;
    const double a121 = 39815761./17514443.;
    const double a124 = -(3457480.+960905.*sq6)/551636.;
    const double a125 = -(3457480.-960905.*sq6)/551636.;
    const double a126 = -844554132./47026969.;
    const double a127 = 8444996352./302158619.;
    const double a128 = -2509602342./877790785.;
    const double a129 = -28888795297996250./3199510091356783.;
    const double a1210 = 226716250./18341897.;
    const double a1211 = 1371316744./2131383595.;
    //@}

    /*! @name Weights (output) */
    //@{
    const double b1 = 104257./1920240.;
    const double b6 = 3399327./763840.;
    const double b7 = 66578432./35198415.;
    const double b8 = -1674902723./288716400.;
    const double b9 = 54980371265625./176692375811392.;
    const double b10 = -734375./4826304;
    const double b11 =  171414593./851261400.;
    const double b12 = 137909./3084480.;
    //@}

    /*! @name Weights (err3) */
    //@{
    const double c1 = 116092271./8848465920.;
    const double c6 = -1871647./1527680.;
    const double c7 = -69799717./140793660.;
    const double c8 = 1230164450203./739113984000.;
    const double c9 = -1980813971228885./5654156025964544.;
    const double c10 = 464500805./1389975552.;
    const double c11 = 1606764981773./19613062656000.;
    const double c12 = -137909./6168960.;
    //@}

    /*! @name Weights (err5) */
    //@{
    const double d1 = -364463./1920240.;
    const double d6 = 3399327./763840.;
    const double d7 = 66578432./35198415.;
    const double d8 = -1674902723./288716400.;
    const double d9 = -74684743568175./176692375811392.;
    const double d10 = -734375./4826304.;
    const double d11 = 171414593./851261400.;
    const double d12 = 69869./3084480.;
    //@}
    //@}


};

#endif // DYN_ODE_DOPRI853_H
